MATHEMATICAL MODELING OF PATTERN FORMATION IN SUB- 
AND SUPPERDIFFUSIVE REACTION-DIFFUSION SYSTEMS 
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Abstract. This paper is concerned with analysis of coupled fractional reaction-diffusion equa- 
tions. It provides analytical comparison for the fractional and regular reaction-diffusion systems. As 
an example, the reaction-diffusion model with cubic nonlinearity and Brusselator model are consid- 
ered. The detailed linear stability analysis of the system with cubic nonlinearity is provided. It is 
shown that by combining the fractional derivatives index with the ratio of characteristic times, it is 
possible to find the marginal value of the index where the oscillatory instability arises. Computer 
simulation and analytical methods are used to analyze possible solutions for a linearized system. 
A computer simulation of the corresponding nonlinear fractional ordinary differential equations is 
presented. It is shown that the increase of the fractional derivative index leads to the periodic solu- 
tions which become stochastic at the index approaching the value of 2. It is established by computer 
simulation that there exists a set of stable spatio-temporal structures of the one-dimensional system 
under the Neumann and periodic boundary condition. The characteristic features of these solutions 
consist in the transformation of the steady state dissipative structures to homogeneous oscillations 
or space temporary structures at a certain value of fractional index. 
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1. Introduction. Reaction-diffusion (RD) systems are inherent in many branches 
of physics, chemistry, biology, ecology etc. The review of the theory and applications 
of reaction-diffusion systems one can find in books and numerous articles (See, for 
example [H [l[3l[l[5lil[7l[8l[9l[T0l[n]). The popularity of the RD system is driven 
by the underlying richness of the nonlinear phenomena, which include stationary and 
spatio-temporary dissipative pattern formation, oscillations, different types of chemi- 
cal waves, excitability, bistability etc. The mechanism of the formation of such type 
of nonlinear phenomena and the conditions of their emergence have been extensively 
studied during the last couple decades. Although the mathematical theory of such 
type of phenomena has not been developed yet due to the essential nonlinearity of 
these systems, from the viewpoint of the applied and experimental mathematics, the 
pattern of possible phenomena in RD system is more or less understandable. 

In the recent years, there has been a great deal of interest in fractional reaction- 
diffusion (FRD) systems [Il[Il[Il[ini[Tll|I71[IllIl[ini[2I] which from one side 
exhibit self-organization phenomena and from the other side introduce a new param- 
eter to these systems, which is a fractional derivative index, and it gives a greater 
degree of freedom for diversity of self-organization phenomena. At the same time, the 
process of analyzing such FRD systems is much complicated from the analytical and 
numerical point of view. 

In this article, we consider two coupled reaction-diffusion systems: 

The first one is the classical system 
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(1.1) ^^ dn,{x t) ^ ;2^2^^(^^^) _ W{n,,n2,A), 

ot 

(1.2) ^^ dn2ix,t) ^ ^2^2^^(^^^) _ 

where x,t € — ■^■,ni{x,t), 712(0;, <) e M - two variables, W ,Q € R are the 

nonlinear sources of the system modeling their production rates, ti,T2,1, L,G M. - 
characteristic times and lengths of the system, ^ S M - is an external parameter 
And the other model is the fractional RD system 

(1.3) ^ ;2^2^^(^^^) _ w{n,,n2,A), 

(1.4) r2^^^^^l^=L'V^n2{x,t)-Qim,n2,A) 

with the same parameters and fractional derivatives ^ g^f on the left hand side 
of equations (|1.3p . (|1.4p instead of standard time derivatives, which are the Caputo 
fractional derivatives in time of the order < a < 2 and are represented as [531 IM] 

d" , , 1 } n(™)(T) 

ai" ^ ' r(m ~a) J (t - r)"+i-™ 


The article is devoted to the second problem and the first one we need for com- 
paring the obtained results with classical one. 

Equations p.3|) . (|1.4p at a = 1 correspond to standard RD system described by 
equations (|l.ip . (|1.2p . At a < 1, they describe anomalous sub-diffusion and at a > 1 - 
anomalous superdiffusion . 

In this paper, we always assume that the following conditions are fulfilled on the 
boundaries 0;lx'- 

(i) Neumann: 

(1.5) dni/dx\x=o = drii / dx\3:^i^ = 

(ii) Periodic: 

(1.6) ni{t,0) = n.i{t,lx), dni/dx\.j,^o ^ dni/dx\x=i^, 
where i = 1,2. 

2. Linear stability analysis. Stability of the steady-state constant solutions 
of the system (|1.3p . p.4p correspond to homogeneous equilibrium states 

(2.1) Wini,n2,A) = 0, Q{ni,n2,A) = 

can be analyzed by linearization of the system nearby this solution. In this case the 
system (|1.3p()1.4p can be transformed to hnear system 

d°'u(x, t) T , s 

(2.2) ' ^Au(x,t), 
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where \i{x,t) = a, ) ,( , ^ = / rr 2^2 \ , an = 

^ ' ' \ An2ix,t) y' ^ -a2i/r2 (L^V-^ - a22)r2 

W^^, ai2 = VF^2' '^21 = Q^n, a22 = Qn2 (^11 derivatives are taken at homoge- 
neous equihbrium states (condition (|2.ip ). By substituting the solution u{x,t) = 

An^(t) ) coskx,k — — 1,2... into FRD system (ll.3p . (|1.4p we can get the 

system of hnear ordinary differential equations (|2.2p with the matrix A determined 
by the operator A, the stabihty conditions of which are given by eigenvalues of this 
matrix. 

Let us analyze the stability of the solution (|2.1[) of the linear system with integer 
derivatives and find the conditions of this instabihty (See, for example: [H |21 SI IS])- 
We repeat this process in order to compare the results obtained with the results of 
the fractional RD system considered in this article. In this case, by searching for 

the solution of the linear system in the form \\.{x,t) = ^ ^ exp(At) cos fcx, we 

get a homogeneous system of linear algebraic equations for constants Ani, An2. The 
solubility of this system leads to the characteristic equation 



(2.3) det |A/ - A| = 0, 

where 



(2.4) A 



{Pk"^ + aii)/Ti axijTx 

021/7-2 (i^fc^ + a22)/T2 



/ is the identity matrix. As a result, the characteristic equation takes on a form of 
a simple quadratic equation — trA\ + det A = 0. 

The linear boundary value problem for RD system , (II. 2p is unstable according 
to inhomogeneous wave vectors /c 7^ if (irA < 0, dot A < 0) 

(2.5) an < -[{1/Lfa22 + 2{l/L){a22aii - 012021)^/'], 



^2(011 + fc^/^) + ri(a22 + fc^L^) > 0, 022011 - 012021 > 

(Turing Bifurcation) and according to homogenous (fc = 0) fluctuations (Hopf Bifur- 
cation) if {trA > 0, det ^ > 0) 

(2.6) r2aii -t- Tia22 < 0, 022O11 - 012O21 > 0. 

For analyzing the equations (|1.3p , p.4p let us also consider a linear system obtained 
near the homogeneous state (12. ip . As a result, the simple linear transformation can 
convert this linear system to a diagonal form 



where C is a diagonal matrix of A: C = P^^AP = (^^q ) ' ^^^^^^^^^^^ ^^^"^ 
are determined by the same characteristic equation (|2.3p with matrix (|2.4p . Ai,2 = 
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(a) (b) 



Fig. 2.1. Schematic view of the marginal curve describing fixed points for two-dimensional 
vector field - (a) the marginal value of a - (b). 



^{trA±y/tr^A - 4det^), r]{t) = P^^ An (ij j ' ^ matrix of eigenvectors 

of matrix A. 

In this case, the solution of the vector equation (|2.7p is given by Mittag-Leffler 
functions [23l EH [23 ESI [22] 



(2.8) An,{t) ^ V / An,(0) = £;^(A,t")An,(0), z = l,2. 

Using the result obtained in the papers [TSl [35] , we can conclude that if for any of the 
roots 

(2.9) \arg{Xi)\ < aTr/2 

the solution has an increasing function component then the system is asymptotically 
unstable. 

Analyzing the roots of the characteristic equations, we can see that at 4det A — 
tr^A > eignevalues are complex and can be represented as 

Ai^2 = ^(.^^^ ± iV4det A - tr'^A) = Xjie ± iXim- 

The roots A1.2 are complex inside the parabola fFigure ICT a)) and the fixed points 
are the spiral source {trA > 0) or spiral sinks {trA < 0). The plot of the marginal 
value a : a — aQ — ^\arg{Xi)\ which follows from the conditions (j2.9p is given by the 
formula 



(2.10) ao 



_ ( f arctan ^4det A/tr2^- 1, trA > 0, 
" \ 2 - I arctan 4 det A / tr'^ A - 1, trA < 
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and is presented in the Figure I2.ir a) below the parabola in the coordinate system 
{trA, a). 

Let us analyze the system solution with the help of the Figure [Ola) . Consider 
the parameters which keep the system inside the parabola. It is a well-known fact, 
that at a = 1 the domain on the righthand side of the parabola {trA > 0) is unstable 
with the existing limit circle, while the domain on the left hand side {trA < 0) is 
stable. By crossing the axis trA — the Hopf bifurcation conditions become true. 
In the general case of a : < a < 2 for every point inside the parabola there 
exists a marginal value of ao where the system changes its stability. The value of a 
is a certain bifurcation parameter which switches the stable and unstable state of the 
system. At lower a : a < = ^\arg{Xi)\ the system has oscillatory modes but they 
are stable. Increasing the value of a > ao = ;||ar(7(Ai)| leads to instability. As a 
result, the domain below the curve ap, as a function of trA is stable and the domain 
above the curve is unstable. 

The plot of the roots, describing the mechanism of the system instability, can be 
understood from the Figure [2TT b) where the case ao > 1 is described. In fact, having 
complex number Ai with ReXi < at a —> 2 it is always possible to satisfy the condi- 
tion |arg(Ai)| < Q;7r/2, and the system becomes unstable according to homogeneous 
oscillations fFigure lCT b)). The smaller is the value of trA, the easier it is to fulfill 
the instability conditions. 

In contrast to this case, a complex values of Ai, with ReXi > lead to the system 
instability for regular system with a = 1. However fractional derivatives with a < 1 
will stabilize the system if a < ao = ^\arg{Xi)\. This makes it possible to conclude 
that fractional derivative equations with a < 1 are more stable that their integer 
twinges. 

3. Solutions of the coupled fractional ordinary differential equations 
(FODEs). Our particular interest here is the analysis of the specific non-linear sys- 
tem of FRD equations. We consider two very well-known examples. The first one is 
the RD system with cubical nonlinearity [21 [H [7] which probably is the simplest one 
used in RD systems modeling 

(3.1) W{ni,n2) = -ni+nl/3 + n2, Q = n2 - /3ni - A. 

The second example is known as Brusselator model [1] and it describes certain chem- 
ical reaction-diffusion processes with a pair of variables whose concentrations are 
controlled by nonlinearities 

(3.2) W{ni,n2) ^ -A+ {f3 + l)ni - nln2, Q{ni,n2) = -(ini + n\n2, 

Let us first consider the coupled fractional ordinary differential equations (FODEs) 
with nonlinearities (|3.ip and analyze the stability conditions for such systems. The 
plot of isoclines for the system (|3.ip is represented on Figure ISTlT a). In this case, 
for homogeneous solution, which can be determined from the system of equations 
W = Q = Q, IS the solution of cubic algebraic equation {[3 — l)ni +ri\/i + A — Q. 
Simple calculation makes it possible to write the expressions required for analysis A = 



It is easy to see that if the value of ti/t2, in certain cases, is smaller than 1, the 
instability conditions ( trA > 0) lead to Hopf bifurcation for regular system (a = 1) 
[U [H [21 [31 [S] . In this case, the plot of the domain, where instability exists, is shown 
on the Figure [XTT bV 
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The linear analysis of the system for a — I shows that, if ti/t2 > 1, the solution 
corresponds to the intersections of two isoclines, and it is stable. The marginal curve, 
separating stability and instability domains, is given on Figure [3?lT b). The smaller 
is the ratio of ri/T2, the wider is the instability region. Formally, at ti/t2 — > 0, 
the instability region in ni coinsides with the interval (—1, 1) where the null isocline 
W{ni,n2) = has its increasing part. The maximum value of the curve Ti/T2(ni) 
corresponds to the value ti/t2 — 1 where the system is neutrally stable. These results 
are very widely known in the theory of nonlinear dynamical systems [I] [31 [31 SI [S] ■ 

In the FODEs the conditions of the instability change l|2.9p . and we have to 
analyze the real and the imaginary part of the existing complex eigenvalues, especially 
the equation: Adet A- tr'^ A = 4((/3 - 1) + n?)/TiT2 - ((1 - nj)/ri - l/r2)^ > 0. In 
fact, with the complex eigenvalues, it is possible to find out the corresponding value 
of a where the condition (|2.9p is true. We show that this interval is not correlated 
with the increasing part of the null isocline of the system. Indeed, omitting simple 

calculation, we can write an equation for marginal values of ni : rif — 2{l + ^)nf + ^ — 

2^(2/3— 1) + 1 = 0, and expression nf = 1 + ^ ±2^/3^ estimates the maximum and 

minimum values of ni where the system can be unstable at certain value of a = ag. 
For example, examine the domain of the FODEs where the eigenvalues are complex 
for fixed value Ti/r2, for example, consider ri = T2 = 1 and /3 = 2. In this case, 
trA = -nl,detA = 1 + nf,4det A - tr'^A = 4 + inj -nj>0, which immediately 

leads to the region of existing of complex roots — V 2 + 2y/2 <ni< V2 + 2V2 and we 
can conclude that the instability region, due to the fractional order of the derivatives, 
can be much wider than the same region for (a ~ 1). The dependance of the value of 
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a on the value (|2.10p is given on Figure [STlT c). 

In this case the plot is obtained at ti/t2 — 1 and that is why the marginal 
instability curve is determined for a > 1 (for a < 1 the system is stable). On this 
figure the domain below the curve corresponds to stability and above it to instability 
conditions. 

Similar analysis can be provided for trA > where (ti/t2) is smaller than unity. 
In this case, the instability conditions are also not correlated with the increasing part 
of the null isocline Wi(ni,n2). In this case the plot of a will be start not from 1 but 
from certain value smaller than unity. However, it is much better for understanding 
to get marginal curve ao as a function of ti/t2. 

Let us analyze the dependence of ao on the parameter ti/t2 where the system 
changes its stability. As it is easy to see from Figure [3TT a).(b) that the easiest way 
to reach instability domain is realized at ^ = when two isoclines are intersected 
themselves at the point (0,0) and this corresponds to maximum of the curve on 
Figure I3.1f b) . Let us consider the dependence of this marginal value of ao on the 
parameter ti /t2 at ^ = and determine the plot of this maximum as a functions 
of ao on ti/t2. Such curve for the given model is represented in the Figure ISTT d). 
This curve obtained from (|2.10p corresponds to the dependence of a on trA (Figure 
I2.1f a)). Below this curve the system is unstable, and above it - it is stable. We 
may therefore focus our attention on the general form of this curve. At sufficiently 
small value ti/t2 oscillatory instability is valid even for small a < 1. In contrast, at 
a > 1, the instability conditions could have place even for those cases when ti/t2 is 
sufficiently large. This means that fractional differential equations, by corresponding 
combination of the parameter ti /t2 , can be stable or unstable practically in all the 
region 1 < a < 2. 

It should be noted that, even if the eigenvalues are not complex (A/m — 0), the 
systems with fractional derivatives can poses oscillatory damping oscillations. Such 
situation takes place when 4 det A — tr^A < 0, trA < 0, det A > and two eigenvalues 
are real and less than zero. In this case, at 1 < a < 2 steady state solutions of the 
system are stable and any perturbations are damping. Such system was considered, 
for example, in the article |31| . where an analytical solution for fractional oscillator is 
obtained. Namely with this case we start analyzing possible solutions. 

Several examples of linear FDOEs were solved analytically by Adomnian decom- 
position method, as well as numerically. The obtained solutions were compared with 
the analytical solutions obtained by Mittag-Leffier functions ()2.8|) and the results of 
the work [31] . 

Three different solutions are plotted on the same Figure 13.21 The results of the 
computer simulation of linear ordinary differential equations for stable system - a) 
and unstable system - b) of one variable ni are given on the time domain interval 
[0,30]. We can see that in the stable domain a < ao the oscillations are damping, 
and at a > ao they grow exponentially. So, the ordinary differential equations (space 
derivatives are equal to zero) of the system have two different modes. The solution is 
either asymptotically stable or unstable. 

It is a statement, that FODEs are at least as stable as their integer order coun- 
terparts [24]. At the same time, we have established here that the dynamics of the 
FODEs can be much more complicated than that of the equations with integer order 
[551 [3H [Mj . Our task here is to clarify these statements by finding out not only the 
conditions of the bifurcation but also the real time dynamics of FODEs. 

The developed technique of Adomian decomposition [571 [Ml [30] is a powerful 
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method for solving the system of ordinary differential equations. The effectiveness 
of this method was demonstrated for solving linear fractional differential equations, 
and we used this method here to test the computer program and to compare these 
analytical results with Mittag-Leffler functions. 

First, we consider the case when the system is always asymptotically stable at 
a < ao = 2, for example a = 1.7 (Figure [3211a)). Three solutions obtained by these 
different methods practically do not differ from each other and they show oscillatory 
damping time behavior (the model described in [31]). Analytical results obtained by 
Adomian decomposition is given by the formula 



which is plot on the Figure I3.2f a) by empty circled curve (we truncate the series at 
k=60). At small time interval (for example t < 7) this series can be represented by 
next expansion 

ni = 1 - 0.6473808267 • i^'^ + 0.9865725648 • 10"^ • t^-^ - 0.7019911214 • 10"^ • t^'^ 
+0.296127716 • 10"^ • t^-^ - 0.838275933 • 10"^ • t^-^ + 0.171848173 • 10"^ • t^°-'^ 
-0.2686528893 • 10"^ • i"-^ + 0.3324725330 • 10"^° • t"'^ 
-0.3350625811 • 10"^^ • t^^-^ + 0.2811457252 • 10"^"* • . 

Despite the considered system, in this particular case, cannot be unstable (Figure 
I3.2r a')') because the values of A are real and negative \arg{Xi)\ = tt > aTr/2, in the 
limit at a = 2 we have a regular linear oscillator [31j , analytical solution of which we 
get from Adomian decomposition in the form of power series. The terms of this series 
completely coincides with the expansion of linear oscillator solution cos{x). 

It should be noted that all damped oscillations solutions of the linear oscillator 
of fractional order have similar plots. 

Quite different is the dynamics of FODEs for a > uq when oscillatory mode 
becomes unstable and leads to increasing oscillations at a > ao (Figure [32Kb)). In 
this case, the analytical solution obtained by Adomian decomposition for a — 1.3 
looks like 

m = 0.67 - 0.1716397314 • t^-^ - 0.05663459553 • t^-^ + 0.008906834532 • t^-^ 

+0.00006991997714 • t^-^ - 0.00004692163249 • t^-^ + 0.1292316598 • 10"^ • f-^ 
+0.5298004324 • 10"^ • t^-^ - 0.2781417477 • 10"^ • t^°-* 
+0.3895105489 • 10~" • t"'^ + 0.1811701048 • 10"" • + ... 

The result of taking into account 50 terms in Adomian decomposition expansion 
makes it possible to represent the solution in the time interval till t = 30 is presented 
on the Figure [n^^b). As a matter of fact, such representation has rather theoretical 
sense because our system is essentially nonlinear and does not allow this type of solu- 
tions. For the initial stage dynamics or for the linear system, Adomian decomposition 
method is very effective. The application of Adomian decomposition to the nonlinear 
FDOEs is not successful. Taking into account 10 terms makes it possible to find out 
the solution on time interval t = 6 — 7 and at higher values of t the discrepancy 
between the numerical solution and the one, obtained by Adomian decomposition, 
increases rapidly. 



(3.3) 
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(a) (b) 

Fig. 3.2. The damped oscillator solution for = —u at a = 1.7, u{0) = 1 - (a), and 
the increasing time domain oscillations of the linearized system il.3{l , [L4\ l for ni and a = 1.3, .4 = 

-0.1,13 = l,Ti = T2 = l,ni(0) = {-3y4)i,n2(0) = {-3A)^ + A - (b). Small filled circled line 
- Mittag- Leffier solution, empty circled line - Adomnian decomposition solution, solid gray line - 
numerical solution 



To demonstrate the nontrivial properties of FDOEs, here we consider the nonhn- 
ear dynamics of the above mentioned nonUnear fractional differential equations. 

For a > ao small perturbation of the steady state solution, due to the memory 
inherent in fractional derivatives survive in the process of evolution and grow in ampli- 
tude while nonlinear terms of the system p.ip restrict the value of these oscillations. 
In this case, time dependence of the variables corresponds to the oscillatory solution. 
(The phase portrait and isochnes are presented on Figure [2!3Ia)-(e)). 

Analyzing the phase trajectory of the FODEs, we can see that the amplitude of the 
oscillations increases with increasing a. At a approaching 2, the oscillations become 
more complicated and at a = 2 they look more quasichaotic. The time dependence 
of this solution is represented in Figure ISTST f). 

Brusselator system with nonlinearities (|3.2p has quite similar behavior. Calculat- 
ing, for example, at ^ = 2, /? = 2, ri = r2 = 1 the value of ao we find that uq — 1.54. 

The phase portrait and isoclines for ()3.2p are presented in Figure [SH] (a)- (c). At 
a < 1.5, we obtain steady state solution and at a > 1.5 - the steady state oscillation. 
The increase of the value a leads to the complication of the phase paths, and the two- 
dimensional phase portrait looks much more complicated (Figure [3T4Ub)-(c)). The 
attractor of the system of the two coupled nonlinear differential equations gets the 
features of strange attractor and at a ^ 2 it corresponds to the attractor of the fourth 
order differential equations determined by the nonlinearities p.2p . 



4. Computer simulation of pattern formation. This section contains a dis- 
cussion of the results of the numerical study of the initial value problem of the system 
(|1.3p(|1.4p . The system with corresponding initial and boundary conditions was inte- 
grated numerically using the explicit and implicit schemes with respect to time and 
centered difference approximation for spatial derivatives. The fractional derivatives 
were approximated using two different schemes on the basis Riemann-Liouville defi- 
nition : Ll-scheme for < a < 1, L2-scheme for 1 < a < 2 (see below and [35]), as 
well as the scheme on the basis of Grunwald-Letnikov definition for < a < 1 and 
1 < a < 2 24J . In other words, for the system of n fractional RD equations 



(4-1) ' £j^ + /j("i, j = l,n, 
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(e) (f) 

Fig. 3.3. Two dimensional phase portrait (a)-(e) and time domain oscillations corresponding 
to plot (e) - (f) of the system 11. 3\ with nonlinearities 13. l\l for A = — 0.1,/3 = 1,ti = T2 = 

l,i = L = 0. (a)-a=1..3, (b) - a=1.7, (c) - a=1.8, (d) - a=1.90, (e) - a=2.0, time domain 
oscillations (f) - a=2.00 



where Tj , dj , fj - certain parameters and nonlinearities of the RD system correspond- 
ingly, we used the next numerical schemes: 
Ll-scheme 



K,,-i-24. + "7,^+i)-/jK.' 
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(c) 



(d) 



Fig. 3.4. Two dimensional phase portrait (a)-(c) and time domain oscillations corresponding 
to plot (c) - (d) of the system 11. 3\ ), fj.^P with nonlinearities IS.SXl for A. = 2, l3 = 2, ti = T2 = 
1,/ = L = 0. (a)-a=1.5, (b)-a=1.6, (c) - a=1.7, time domain oscillations (f) - a=1.7 



Tj-(A^)-"^ (a,) _ 1 
r(2-a,) 



k'-"' - (fc - 1) 



l—Oii 



L2-scheme 



r(p - + 1) 9if 



;=2 p=o 



_ r,-(A^)-"^ 
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= - 2fc^-"^ + 3(fc - 1)2-"^- - (fc - 2)2-"^ , 



and G-L scheme 



r,(Ax) 



c^) = i, c^) = c|!i)(^i-— ^j, ; = i,2,... 

where = Uj{xi,tk) = Uj{iAx, kAt), m= [a]. 

The appUed numerical schemes are implicit, and for each time layer they are 
presented as the system of algebraic equations solved by Newton- Raphson technique. 
Such approach makes it possible to get the system of equations with band Jacobian for 
each node and to use the sweep method for the solution of linear algebraic equations. 
Calculating the values of the spatial derivatives and corresponding nonlinear terms 
on the previous layer, we obtained explicit schemes for integration. Despite the fact 
that these algorithms are quite simple, they are very sensitive and require small steps 
of integration, and they often do not allow to find numerical results. In contrast, 
the implicit schemes, in certain sense, are similar to the implicit Euler's method, and 
they have shown very good behavior at the modeling of fractional reaction-diffusion 
systems for different step size of integration, as well as for nonlinear function and the 
power function of fractional index. Moreover, by modeling according to this algorithm 
system (ll.ip . (|1.2p . we have observed that these results fully match the results obtained 
prior. 

It should be noted that the definition of the fractional derivative in Grunwald- 
Letnikow form is equivalent to the one in Ricmann-Liouvillc method, but for numerical 
calculations it is much more flexible. 

We have considered here the kinetics of formation of dissipative structures for 
different values of a. These results are presented on Figures W7\\ and W?]\ 

The simulations were carried out for a one-dimensional system on an equidistant 
grid with spatial step h changing from = 0.1 to 0.01 and time step At changing from 
0.001 to 0.1. We used imposed Neuman (jl.Sp or periodic boundary conditions p.6p . 
As the initial condition, we used the uniform state which was superposed with a small 
spatially inhomogeneous perturbation. 

The systems have rich dynamics, including steady state dissipative structures, 
homogeneous and nonhomogeneous oscillations, and spatiotemporal patterns. In this 
paper, we focus mainly on the study of general properties of the solutions depending 
on the value of a. 

As discussed in Section 2, there are two different regions in parameter A, where 
the system can be stable or unstable. In the case of a = 1 the steady state solutions 
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Fig. 4.1. Numerical solution of the fractional reaction- diffusion equations il.3\ ), (l.4\ l with 
nonlinearities iS. IV . Dynamics of variable ni on the time interval (0,50) for Ix = ii, A = —0.1,13 = 
l,Ti = T2 = 1, «2 = 0.05, = 1; (a) - a=0.1, (b) - a=0.5, (c) - a=0.99, (d) - a=1.5, (e) - a= 
1.6, (f) ~a=1.8 

in the form of nonhomogcneous dissipative structures are inherent to unstable region 
Til G (—1, !)■ Figures [331 a-)- (c) show the steady state dissipative structure formation 
and Figures [17T] fd)-(f) present the spatio-temporal evolution of dissipative structures, 
which eventually leads to homogeneous oscillations. 

On the Figure HTlT aVfd). the value a increases from 0.1 to 1.5 and on this whole 
interval the structures are in steady state. This is due to the case a < ao, the 
oscillatory perturbations are damping, and we can see that small oscillations are at 
the transition period ni. With increasing a, the steady state structures change to the 
spatio-temporary behavior (Figure iLlf eVfl)). 

The emergence of homogeneous oscillations, which destroy pattern formation 
(Figure [131e),(f)) has deep physical meaning. The matter is that the stationary 
dissipative structures consist of smooth and sharp regions of variable ni, and the 
smooth shape of n2 . The linear system analysis shows that the homogeneous distribu- 
tion of the variables is unstable according to oscillatory perturbations inside the wide 
interval of ni, which is much wider then interval (—1, 1). At the same time, smooth 
distributions at the maximum and minimum values of ni are ±\/3 correspondingly. 
In the first approximation, these smooth regions of the dissipative structures resem- 
ble homogeneous ones and are located inside the instability regions. As a result, the 
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(c) (f) 

Fig. 4.2. Numerical solution of the fractional reaction- diffusion equations il.3\ ), (l.4f with 
nonlinearities (3^. Dynamics of variable n\ on the time interval (0,20) for Ix = 10, A = 2,/3 = 
2,Ti = T2 = 1, = 0.1, = 10; (a) - a=0.1, (b) - a=0.5, (c) - a=0.99, (d) - a=1.5, (e) - a = 
1.6, (f) -a=1.8 

unstable fluctuations lead to homogeneous oscillations, and the dissipative structures 
destroy themselves. We can conclude that oscillatory modes in such type FODEs 
have a much wider attraction region than the corresponding region of the dissipative 
structures. 

For a wide range of the parameters a, the numerical solutions of the Brusselator 
problem show similar behavior (Figure [42r a)-(d)). The stationary solutions emerge 
practically in the same way. At small a, we see aperiodic formation of the structures, 
and approaching ao, the damping oscillations of the dissipative structures arise. At 
ao = 1-7 certain non-stationary structures arise (Figure [12{e)). In this case, the 
dissipative structures look quite similar to those we obtained for regular system [36l 
l37] . The increase of a leads to a larger amplitude of pulsation. All these patterns are 
robust with respect to small initial perturbations. The further increase of a leads to 
spatially temporary chaos (Figure |42tf)). 

In the contrast to previous case, such nonhomogeneous behavior is stable and 
does not lead to homogeneous oscillations. The matter is that in Brusselator model, 
the dissipative structures are much grater in amplitude and do not have smooth 
distribution at the top. 

It should be noted that the pulsation phenomena of the dissipative structures is 
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closely related to the oscillation solutions of the ODE (Figures 13.31 13. 4p . Moreover, 
the fractional derivative of the first variable has the most impact on the oscillations 
emergence. It can be obtained by performing a simulation where the first variable 
is a fractional derivative and the second one is an integer. It should be emphasized 
that the distribution of n2, within the solution, only shows a small deviation from the 
stationary value (that is why this variable is not represented in the figures). 

5. Conclusion. In this article we developed a linear theory of instability of 
reaction diffusion system with fractional derivatives. The introduced new parameter 
- marginal value ao plays the role of bifurcation parameter. If the fractional derivative 
index a is smaller than ap, the system of FODEs is stable and has oscillatory damping 
solutions. At a > OfQ, the FODEs becomes unstable, and we obtain oscillatory or even 
more complex - quasi chaotic solutions. In addition, the stable and unstable domains 
of the system were investigated. 

By the computer simulation of the fractional reaction-diffusion systems we pro- 
vided evidence that pattern formation in the fractional case, at a less than a certain 
value, is practicably the same as in the regular case scenario a = 1. At a > aoj the 
kinetics of formation becomes oscillatory. At a = aoj the oscillatory mode arises and 
can lead to homogeneous or nonhomogeneous oscillations. In the last case scenario, 
depending on the parameters of the medium, we can see a rich variety of spatiotem- 
poral behavior. 
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